function [lb q] = compute_expected( tau , k )

tau = tau( : , 1:k );
digamma_tau = psi( tau );
digamma_sum = psi( sum( tau ) );

digamma_tau1_cumsum = [ 0 cumsum( digamma_tau( 1 , : ) ) ] ;
digamma_sum_cumsum = cumsum( digamma_sum );
tmp = digamma_tau( 2 , : ) + digamma_tau1_cumsum(1:k) - digamma_sum_cumsum;
q = exp( tmp - max(tmp) );
q = q / sum(q);


lb = sum( q .* ( tmp - log( q ) ) );

return
